home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlasq2.f < prev    next >
Text File  |  1996-07-19  |  9KB  |  269 lines
  1.       SUBROUTINE DLASQ2( M, Q, E, QQ, EE, EPS, TOL2, SMALL2, SUP, KEND,
  2.      $                   INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            INFO, KEND, M
  11.       DOUBLE PRECISION   EPS, SMALL2, SUP, TOL2
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   E( * ), EE( * ), Q( * ), QQ( * )
  15. *     ..
  16. *
  17. *     Purpose
  18. *     =======
  19. *
  20. *     DLASQ2 computes the singular values of a real N-by-N unreduced
  21. *     bidiagonal matrix with squared diagonal elements in Q and
  22. *     squared off-diagonal elements in E. The singular values are
  23. *     computed to relative accuracy TOL, barring over/underflow or
  24. *     denormalization.
  25. *
  26. *     Arguments
  27. *     =========
  28. *
  29. *  M       (input) INTEGER
  30. *          The number of rows and columns in the matrix. M >= 0.
  31. *
  32. *  Q       (output) DOUBLE PRECISION array, dimension (M)
  33. *          On normal exit, contains the squared singular values.
  34. *
  35. *  E       (workspace) DOUBLE PRECISION array, dimension (M)
  36. *
  37. *  QQ      (input/output) DOUBLE PRECISION array, dimension (M)
  38. *          On entry, QQ contains the squared diagonal elements of the
  39. *          bidiagonal matrix whose SVD is desired.
  40. *          On exit, QQ is overwritten.
  41. *
  42. *  EE      (input/output) DOUBLE PRECISION array, dimension (M)
  43. *          On entry, EE(1:N-1) contains the squared off-diagonal
  44. *          elements of the bidiagonal matrix whose SVD is desired.
  45. *          On exit, EE is overwritten.
  46. *
  47. *  EPS     (input) DOUBLE PRECISION
  48. *          Machine epsilon.
  49. *
  50. *  TOL2    (input) DOUBLE PRECISION
  51. *          Desired relative accuracy of computed eigenvalues
  52. *          as defined in DLASQ1.
  53. *
  54. *  SMALL2  (input) DOUBLE PRECISION
  55. *          A threshold value as defined in DLASQ1.
  56. *
  57. *  SUP     (input/output) DOUBLE PRECISION
  58. *          Upper bound for the smallest eigenvalue.
  59. *
  60. *  KEND    (input/output) INTEGER
  61. *          Index where minimum d occurs.
  62. *
  63. *  INFO    (output) INTEGER
  64. *          = 0:  successful exit
  65. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  66. *          > 0:  if INFO = i, the algorithm did not converge;  i
  67. *                specifies how many superdiagonals did not converge.
  68. *
  69. *  =====================================================================
  70. *
  71. *     .. Parameters ..
  72.       DOUBLE PRECISION   ZERO
  73.       PARAMETER          ( ZERO = 0.0D+0 )
  74.       DOUBLE PRECISION   FOUR, HALF
  75.       PARAMETER          ( FOUR = 4.0D+0, HALF = 0.5D+0 )
  76. *     ..
  77. *     .. Local Scalars ..
  78.       INTEGER            ICONV, IPHASE, ISP, N, OFF, OFF1
  79.       DOUBLE PRECISION   QEMAX, SIGMA, XINF, XX, YY
  80. *     ..
  81. *     .. External Subroutines ..
  82.       EXTERNAL           DLASQ3
  83. *     ..
  84. *     .. Intrinsic Functions ..
  85.       INTRINSIC          MAX, MIN, NINT, SQRT
  86. *     ..
  87. *     .. Executable Statements ..
  88.       N = M
  89. *
  90. *     Set the default maximum number of iterations
  91. *
  92.       OFF = 0
  93.       OFF1 = OFF + 1
  94.       SIGMA = ZERO
  95.       XINF = ZERO
  96.       ICONV = 0
  97.       IPHASE = 2
  98. *
  99. *     Try deflation at the bottom
  100. *
  101. *     1x1 deflation
  102. *
  103.    10 CONTINUE
  104.       IF( N.LE.2 )
  105.      $   GO TO 20
  106.       IF( EE( N-1 ).LE.MAX( QQ( N ), XINF, SMALL2 )*TOL2 ) THEN
  107.          Q( N ) = QQ( N )
  108.          N = N - 1
  109.          IF( KEND.GT.N )
  110.      $      KEND = N
  111.          SUP = MIN( QQ( N ), QQ( N-1 ) )
  112.          GO TO 10
  113.       END IF
  114. *
  115. *     2x2 deflation
  116. *
  117.       IF( EE( N-2 ).LE.MAX( XINF, SMALL2,
  118.      $    ( QQ( N ) / ( QQ( N )+EE( N-1 )+QQ( N-1 ) ) )*QQ( N-1 ) )*
  119.      $    TOL2 ) THEN
  120.          QEMAX = MAX( QQ( N ), QQ( N-1 ), EE( N-1 ) )
  121.          IF( QEMAX.NE.ZERO ) THEN
  122.             IF( QEMAX.EQ.QQ( N-1 ) ) THEN
  123.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  124.      $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
  125.      $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
  126.             ELSE IF( QEMAX.EQ.QQ( N ) ) THEN
  127.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  128.      $              SQRT( ( ( QQ( N-1 )-QQ( N )+EE( N-1 ) ) /
  129.      $              QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) )
  130.             ELSE
  131.                XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX*
  132.      $              SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) /
  133.      $              QEMAX )**2+FOUR*QQ( N-1 ) / QEMAX ) )
  134.             END IF
  135.             YY = ( MAX( QQ( N ), QQ( N-1 ) ) / XX )*
  136.      $           MIN( QQ( N ), QQ( N-1 ) )
  137.          ELSE
  138.             XX = ZERO
  139.             YY = ZERO
  140.          END IF
  141.          Q( N-1 ) = XX
  142.          Q( N ) = YY
  143.          N = N - 2
  144.          IF( KEND.GT.N )
  145.      $      KEND = N
  146.          SUP = QQ( N )
  147.          GO TO 10
  148.       END IF
  149. *
  150.    20 CONTINUE
  151.       IF( N.EQ.0 ) THEN
  152. *
  153. *         The lower branch is finished
  154. *
  155.          IF( OFF.EQ.0 ) THEN
  156. *
  157. *         No upper branch; return to DLASQ1
  158. *
  159.             RETURN
  160.          ELSE
  161. *
  162. *         Going back to upper branch
  163. *
  164.             XINF = ZERO
  165.             IF( EE( OFF ).GT.ZERO ) THEN
  166.                ISP = NINT( EE( OFF ) )
  167.                IPHASE = 1
  168.             ELSE
  169.                ISP = -NINT( EE( OFF ) )
  170.                IPHASE = 2
  171.             END IF
  172.             SIGMA = E( OFF )
  173.             N = OFF - ISP + 1
  174.             OFF1 = ISP
  175.             OFF = OFF1 - 1
  176.             IF( N.LE.2 )
  177.      $         GO TO 20
  178.             IF( IPHASE.EQ.1 ) THEN
  179.                SUP = MIN( Q( N+OFF ), Q( N-1+OFF ), Q( N-2+OFF ) )
  180.             ELSE
  181.                SUP = MIN( QQ( N+OFF ), QQ( N-1+OFF ), QQ( N-2+OFF ) )
  182.             END IF
  183.             KEND = 0
  184.             ICONV = -3
  185.          END IF
  186.       ELSE IF( N.EQ.1 ) THEN
  187. *
  188. *     1x1 Solver
  189. *
  190.          IF( IPHASE.EQ.1 ) THEN
  191.             Q( OFF1 ) = Q( OFF1 ) + SIGMA
  192.          ELSE
  193.             Q( OFF1 ) = QQ( OFF1 ) + SIGMA
  194.          END IF
  195.          N = 0
  196.          GO TO 20
  197. *
  198. *     2x2 Solver
  199. *
  200.       ELSE IF( N.EQ.2 ) THEN
  201.          IF( IPHASE.EQ.2 ) THEN
  202.             QEMAX = MAX( QQ( N+OFF ), QQ( N-1+OFF ), EE( N-1+OFF ) )
  203.             IF( QEMAX.NE.ZERO ) THEN
  204.                IF( QEMAX.EQ.QQ( N-1+OFF ) ) THEN
  205.                   XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+
  206.      $                 QEMAX*SQRT( ( ( QQ( N+OFF )-QQ( N-1+OFF )+EE( N-
  207.      $                 1+OFF ) ) / QEMAX )**2+FOUR*EE( OFF+N-1 ) /
  208.      $                 QEMAX ) )
  209.                ELSE IF( QEMAX.EQ.QQ( N+OFF ) ) THEN
  210.                   XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+
  211.      $                 QEMAX*SQRT( ( ( QQ( N-1+OFF )-QQ( N+OFF )+EE( N-
  212.      $                 1+OFF ) ) / QEMAX )**2+FOUR*EE( N-1+OFF ) /
  213.      $                 QEMAX ) )
  214.                ELSE
  215.                   XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+
  216.      $                 QEMAX*SQRT( ( ( QQ( N+OFF )-QQ( N-1+OFF )+EE( N-
  217.      $                 1+OFF ) ) / QEMAX )**2+FOUR*QQ( N-1+OFF ) /
  218.      $                 QEMAX ) )
  219.                END IF
  220.                YY = ( MAX( QQ( N+OFF ), QQ( N-1+OFF ) ) / XX )*
  221.      $              MIN( QQ( N+OFF ), QQ( N-1+OFF ) )
  222.             ELSE
  223.                XX = ZERO
  224.                YY = ZERO
  225.             END IF
  226.          ELSE
  227.             QEMAX = MAX( Q( N+OFF ), Q( N-1+OFF ), E( N-1+OFF ) )
  228.             IF( QEMAX.NE.ZERO ) THEN
  229.                IF( QEMAX.EQ.Q( N-1+OFF ) ) THEN
  230.                   XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+
  231.      $                 QEMAX*SQRT( ( ( Q( N+OFF )-Q( N-1+OFF )+E( N-1+
  232.      $                 OFF ) ) / QEMAX )**2+FOUR*E( N-1+OFF ) /
  233.      $                 QEMAX ) )
  234.                ELSE IF( QEMAX.EQ.Q( N+OFF ) ) THEN
  235.                   XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+
  236.      $                 QEMAX*SQRT( ( ( Q( N-1+OFF )-Q( N+OFF )+E( N-1+
  237.      $                 OFF ) ) / QEMAX )**2+FOUR*E( N-1+OFF ) /
  238.      $                 QEMAX ) )
  239.                ELSE
  240.                   XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+
  241.      $                 QEMAX*SQRT( ( ( Q( N+OFF )-Q( N-1+OFF )+E( N-1+
  242.      $                 OFF ) ) / QEMAX )**2+FOUR*Q( N-1+OFF ) /
  243.      $                 QEMAX ) )
  244.                END IF
  245.                YY = ( MAX( Q( N+OFF ), Q( N-1+OFF ) ) / XX )*
  246.      $              MIN( Q( N+OFF ), Q( N-1+OFF ) )
  247.             ELSE
  248.                XX = ZERO
  249.                YY = ZERO
  250.             END IF
  251.          END IF
  252.          Q( N-1+OFF ) = SIGMA + XX
  253.          Q( N+OFF ) = YY + SIGMA
  254.          N = 0
  255.          GO TO 20
  256.       END IF
  257.       CALL DLASQ3( N, Q( OFF1 ), E( OFF1 ), QQ( OFF1 ), EE( OFF1 ), SUP,
  258.      $             SIGMA, KEND, OFF, IPHASE, ICONV, EPS, TOL2, SMALL2 )
  259.       IF( SUP.LT.ZERO ) THEN
  260.          INFO = N + OFF
  261.          RETURN
  262.       END IF
  263.       OFF1 = OFF + 1
  264.       GO TO 20
  265. *
  266. *     End of DLASQ2
  267. *
  268.       END
  269.